/* 
 * zx_fft.c 
 * 
 * Implementation of Fast Fourier Transform(FFT) 
 * and reversal Fast Fourier Transform(IFFT) 
 *  
 *  Created on: 2013-8-5 
 *      Author: monkeyzx 
 * 
 * TEST OK 2014.01.14 
 * == 2014.01.14 
 *   Replace @BitReverse(x,x,N,M) by refrence to  
 *   <The Scientist and Engineer's Guide to Digital Signal Processing> 
 */  
  
#include "zx_fft.h"  

  
/* 
 * FFT Algorithm 
 * === Inputs === 
 * x : complex numbers 
 * N : nodes of FFT. @N should be power of 2, that is 2^(*) 
 * === Output === 
 * the @x contains the result of FFT algorithm, so the original data 
 * in @x is destroyed, please store them before using FFT. 
 */  
int fft(TYPE_FFT *x, int N)  
{  
    int i,j,l,k,ip;  
    int M = 0;  
    int le,le2;  
	double sR,sI,tR,tI,uR,uI;  
  
    M = (int)log2(N);  
  
    /* 
     * bit reversal sorting 
     */  
    l = N / 2;  
    j = l;  
    //BitReverse(x,x,N,M);  
    for (i=1; i<=N-2; i++) {  
        if (i < j) {  
            tR = x[j].real;  
            tI = x[j].imag;  
            x[j].real = x[i].real;  
            x[j].imag = x[i].imag;  
            x[i].real = tR;  
            x[i].imag = tI;  
        }  
        k = l;  
        while (k <= j) {  
            j = j - k;  
            k = k / 2;  
        }  
        j = j + k;  
    }  
  
    /* 
     * For Loops 
     */  
    for (l=1; l<=M; l++) {   /* loop for ceil{log2(N)} */  
        le = (int)pow(2,l);  
        le2 = (int)(le / 2);  
        uR = 1;  
        uI = 0;  
        sR = cos(PI / le2);  
        sI = -sin(PI / le2);  
        for (j=1; j<=le2; j++) {   /* loop for each sub DFT */  
            //jm1 = j - 1;  
            for (i=j-1; i<=N-1; i+=le) {  /* loop for each butterfly */  
                ip = i + le2;  
                tR = x[ip].real * uR - x[ip].imag * uI;  
                tI = x[ip].real * uI + x[ip].imag * uR;  
                x[ip].real = x[i].real - tR;  
                x[ip].imag = x[i].imag - tI;  
                x[i].real += tR;  
                x[i].imag += tI;  
            }  /* Next i */  
            tR = uR;  
            uR = tR * sR - uI * sI;  
            uI = tR * sI + uI *sR;  
        } /* Next j */  
    } /* Next l */  
  
    return 0;  
}  
  
/* 
 * Inverse FFT Algorithm 
 * === Inputs === 
 * x : complex numbers 
 * N : nodes of FFT. @N should be power of 2, that is 2^(*) 
 * === Output === 
 * the @x contains the result of FFT algorithm, so the original data 
 * in @x is destroyed, please store them before using FFT. 
 */  
int ifft(TYPE_FFT *x, int N)  
{  
    int k = 0;  
  
    for (k=0; k<=N-1; k++) {  
        x[k].imag = -x[k].imag;  
    }  
  
    fft(x, N);    /* using FFT */  
  
    for (k=0; k<=N-1; k++) {  
        x[k].real = x[k].real / N;  
        x[k].imag = -x[k].imag / N;  
    }  
  
    return 0;  
} 


int zeropadding(float*x,int N1,TYPE_FFT *temp,int N2){
	int i;
	//TYPE_FFT *temp=malloc(N2*sizeof(TYPE_FFT));
	for(i=0;i<N1;i++){
		temp[i].real=x[i];
		temp[i].imag=0;
	}
	for(i=N1;i<N2;i++){
		temp[i].real=0;
		temp[i].imag=0;
	}
	return 0;
	
}
